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We present dynamic measurements of a lattice gas during phase separation, which we use as an 
analogy for self-assembly of equilibrium ordered structures. We use two approaches to quantify the 
degree of 'reversibility' of this process: firstly, we count events in which bonds are made and bro- 
ken; secondly, we use correlation-response measurements and fluctuation-dissipation ratios to probe 
reversibility during different time intervals. We show how correlation and response functions can 
be related directly to microscopic (ir)reversibility and we discuss time-dependence and observable- 
dependence of these measurements, including the role of fast and slow degrees of freedom during 
assembly. 

PACS numbers: 



I. INTRODUCTION 

Self-assembly [JHH] is the spontaneous formation of 
complex ordered equilibrium structures from simpler 
component particles. The range of possible structures 
includes novel crystals @H7]> viral capsids 0|S], and col- 
loidal molecules [T0Tfl"2] . Self-assembly is viewed as a 
promising alternative technology to current fabrication 
techniques, offering a bottom-up approach to design and 
manufacture |13j . For example, experimental work has 
used DNA to generate a backbone for potential nanofab- 
rication [T4"l IT5] ; potential applications of artificial viral 
capsids include inert vaccines and drug delivery |16j ; and 
the potential range of self-assembled structures available 
through control of particle shape and interactions have 
been discussed extensively [2J [51 OH] • Self-assembly has 
even been considered for the potential regeneration of 
human organs and tissues [T7]- However a key challenge 
for design and control of self-assembly processes is that 
even in systems where the equilibrium states are known, 
the conditions which lead to effective assembly remain 
poorly understood and there remains no general theoret- 
ical approach to support experimental advances. 

A major step towards such a theoretical approach has 
been the recognition of the importance of reversibility 
in the assembly process [S} fI5H2I] . As a system evolves 
from an initially disordered state it makes bonds between 
the constituent particles. If structures form which are 
not typical of equilibrium configurations then they must 
anneal before the system arrives at equilibrium. When 
bonds between particles are too strong, thermal fluctu- 
ations are insufficient to break incorrect bonds before 
additional particles aggregate, and the system becomes 
trapped in long-lived kinetically frustrated states. To 
avoid this problem and achieve effective self-assembly, 
bonds must be both made and broken as the system 
evolves in time. In this sense, self-assembling systems 
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are typically reversible on short time scales, even though 
they change macroscopically on long timescales. 

Recent work |19l 12(71 [221 123"] has sought to quantify this 
qualitative idea by making dynamic measurements of re- 
versibility. We will compare two approaches in a simple 
lattice gas model where particles assemble into clusters 
and the effectiveness of self-assembly is identified by the 
presence (or absence) of under-coordinated particles in 
the clusters. Firstly, we count individual bonding and 
unbonding events and compare their relative frequencies 
by defining dynamical flux and traffic observables [23J. 
Then we use measurements of responses to an applied 
field to probe 'reversibility', using fluctuation-dissipation 
theory |24f[2rj] combined with measurements of out-of- 
equilibrium correlation functions. We relate measure- 
ments of fluctuation dissipation ratios (FDRs) |2B] to 
the reversibility of self-assembly, extending the analysis 
of |19l [22] , In particular, we derive a formula for the 
response that elucidates its relation to microscopic re- 
versibility and to measurements of flux and traffic |23| . 
We discuss how these results affect the idea |19l 122] that 
measurements of reversibility on short time scales might 
be used to predict long-time behaviour. 



II. THE MODEL 
A. Definition 

We use an Ising lattice g simple model system 

for self-assembly. While this model appears much sim- 
pler than more detailed models of crystallisation [4, 7 
or viral capsid self-assembly [9j [20], previous work has 
shown that it mirrors many of the physical phenomena 
that occur in model self-assembling systems, particularly 
kinetic trapping [2T1 1251 [27] . 

At low temperatures, the equilibrium state of the lat- 
tice gas consists of large dense clusters of particles, but ef- 
ficient assembly of these clusters requires reversible bond- 
ing in order to avoid kinetic trapping (particularly aggre- 
gation into ramified fractal structures [25]). 
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FIG. 1: a) Exact phase diagram for the lattice gas |29| . At 
high temperature the system equilibrates in a single fluid 
phase while below the binodal it separates into high and low 
density phases, b) Plot of yield 714(f) against reduced temper- 
ature T/eb, for t = 10 4 and 10 6 MCS. For this range of times, 
the yield is maximised at T/eb ~ 0.35. The equilibrium yield 
is shown (labelled 'eqbm'): as t — > 00 the yield approaches 
this result. The binodal at T/eb = 0.547 is shown as a vertical 
line. 



The model consists of a lattice, each site of which may 
contain at most one particle, and the energy of the system 
is 



E 



£b 
2 



(i) 



where the summation is over all particles, £h is the inter- 
action strength, and n p the number of occupied nearest 
neighbours of particle p. We consider N = 1638 parti- 
cles on a square lattice of dimension d = 2 and linear 
size L — 128. The relevant system parameters are the 
density p = N/L d ps 0.1 and the bond strength e b /T (or 
equivalently, the reduced temperature T/eb). The phase 
diagram is shown in Fig. [T] indicating the single fluid 
phase and the region of liquid-gas coexistence which oc- 
curs when sinh 4 (e b /2T) > (1 - (2p - l) 8 )" 1 ^. At the 
density p — 0.1 at which we work, the binodal is located 
at T/eb = 0.547, and below this temperature the system 
separates into high and low density phases at equilib- 
rium (the critical point for the model is at p c = 0.5, 
Tc/tb = 0.567). For the analogy with self-assembly that 
we are considering, we are interested in behaviour at 
fairly low densities, as in (9j H0 1I20H22] , All of the data 
that we show is taken at p — 0.1 but the same qual- 
itative behaviour is found for lower densities too. (At 
higher densities, we find that clusters of particles start 
to percolate through the system. We do not consider 
this regime since it would corresponds to gelation in the 
self-assembling systems, and this is not relevant for the 
systems we have in mind.) 

We use a Monte Carlo (MC) scheme to simulate the 
diffusive motion of particles as they assemble. Our cen- 
tral assumptions are (i) that clusters of n particles diffuse 
with rates proportional to 1/n, and (ii) that bond-making 
is diffusion-limited (that is, bond-making rates depend 
weakly on the bond strength eb/T while bond-breaking 
rates have an Arrhenius dependence on eb/T.) Our sim- 



ulations begin from a random arrangement of particles, 
and we simulate dynamics at fixed bond strength eb/T. 
The dynamics are based on the 'cleaving' algorithm [3D] 
which makes cluster moves in accordance with detailed 
balance and ensures physically realistic diffusion. We 
propose clusters of particles to be moved by picking a 
seed particle at random: each neighbour of that particle 
is added to the cluster with probability 1 — exp (— Aeb/T) 
where A is a parameter that determines the relative like- 
lihood of cluster rearrangement and cluster motion. We 
take A = 0.9 (see below). The cluster to be moved is built 
up by recursively adding neighbours of those particles in 
the cluster: each possible bond is tested once. 

In addition, a maximum cluster size n max is chosen 
for each proposed move, and the move may be accepted 
only if the size of the cluster to be moved is less than 
"■max- We choose n max to be a real number, greater than 
unity, distributed as P(n max > n) — 1/n 2 . The result 
is that clusters of n particles have a diffusion constant 
D cx 1/n, consistent with Brownian dynamics: see |30| 
for full details. 

Having generated a cluster, one of the four lattice di- 
rections is chosen at random, and we attempt to move 
the cluster a single lattice spacing in that direction, re- 
jecting any moves which cause particles to overlap. Fi- 
nally, the proposed cluster move is accepted with prob- 
ability P a . which depends on its energy change AE. 
For AE 7^ 0, we take a Metropolis formula P a = 
min{l,exp[(A - 1)AE/T})} while for AE = we take 
P a = a with a = 0.9 a constant (see below). (Note that 
the definition of the energy in this model means that 
AE/eb is an integer, which ensures that P a is monotonic 
in AE for all values of eb that we consid er: cas es where 
AE/eb is non-integer are discussed in Sec. Ill C| ) Simula- 
tion time, t, is measured in Monte Carlo sweeps (MCS), 
with 1MCS corresponding to N attempted MC moves; 
the time associated with a single attempted MC move is 
therefore St — 1/N. We note that the constants a and 
A within this algorithm may be chosen freely between 
and 1, while still preserving detailed balance. Our MC 
algorithm mimics diffusive cluster motion most closely 
when a and A are close to unity. However, when evalu- 
ating response functions (see below), the MC transition 
rates should be differentiablc functions of any perturbing 
fields: this may not be the case if either a and A is equal 
to unity. We therefore take a — 0.9 and A = 0.9. 



B. Effectiveness of self-assembly 

In studying the lattice gas we have in mind a sys- 
tem self-assembling (or phase separating) into an ordered 
phase such as a crystal, coexisting with a dilute fluid. 
In the ordered phase, particles have a specific local co- 
ordination environment. In the lattice gas model, we 
take this local environment to be that a particle has all 
neighbouring sites occupied. We define the 'assembly 
yield', 714(f), as the proportion of particles with the max- 
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FIG. 2: Snapshots of configurations at three representative 
temperatures. At eb/T = 0.15 (low temperature), long-lived 
fractal clusters survive to form a stable gel-like structure, 
which we identify as kinetically frustrated. At eb/T = 0.5 
(a higher temperature) and for long times, the system fully 
phase separates into a large cluster surrounded by a dilute gas. 
However, the large cluster is not close-packed and a large pro- 
portion of particles remain in small clusters: this is therefore 
poor assembly. At eb/T = 0.35 (intermediate temperature) 
the majority of the particles are in large clusters having few 
defects, which we identify as good assembly. 



imal number of 4 bonds to neighbouring particles. In 
Fig. [T]d we plot the average yield at two fixed times. We 
observe a maximum at a reduced temperature that we 
denote by T* /e\,: the value of T* depends only weakly 
on time t. The system shows three qualitative kinds of 
behaviour. To illustrate this, we show snapshots from 
our simulations in Fig. [2j at three representative bond 
strengths. For strong bonds, there is a decrease in yield 
as the system becomes trapped in kinetically frustrated 
states. For weak bonds, the system is already close to 
equilibrium at t = 10 6 but the final state has a rela- 
tively low value of n±{t). In the language of self-assembly, 
we interpret this as a poor quality product. Based on 
Figs. [T] and [2] we loosely identify the kinetically trapped 
regime as T/eh < 0.2 and the regime of weak bond- 
ing and poor assembly as T/e^ > 0.5. However, we 
emphasise that these regimes are separated by smooth 
crossovers and not sharp transitions. In the following, we 
use T/e b = 0.15,0.35,0.5 as representative state points 
for kinetic trapping, good assembly and poor assembly, 
respectively. 

Thus, despite the simplicity of the lattice gas, it cap- 
tures the kinetic trapping effects and the non-monotonic 
yield observed in more physically realistic model sys- 
tems [9] flQ| Hfj-23 . We emphasise that the kinetically 
trapped states we find are closely related to diffusion- 
limited aggregates |28J, while the 'assembled' states are 
compact clusters. The changes in cluster morphology on 
varying bond strength are discussed further in [27] but in 
this article we concentrate on the dynamical reversibility 
of assembly [18-- 20, 22 and not on the structure of the 



clusters. We include in Fig. [T] the equilibrium behaviour 
for systems of this size, which we obtain by running dy- 
namical simulations starting from a fully phase-separated 
state in which a single large cluster contains all particles. 
As t — > oo, the yield must approach its equilibrium re- 
sult: the laws of thermodynamics state that if we wait 
long enough then the highest quality assembly (or the 
lowest energy final state) will be at the lowest tempera- 
ture. 



III. MEASURING REVERSIBILITY 

A. Flux- Traffic Ratio 

We now turn to measurements of reversibility and 
their relation to kinetic trapping effects and the non- 
monotonic yield shown in Fig. [T| We follow [S3] in con- 
sidering the net rate of energy changing events at a mi- 
croscopic level, the flux, in proportion to the total rate of 
energy changing events, the traffic (see also |3"D - l33"] ). We 
count an event (or 'kink') whenever a particle changes its 
number of neighbours: the number of times that parti- 
cle i increases its number of bonds between times t and 
t + At is K+(t, At), with K~(t, At) the number of times 
the particle decreases its number of bonds. We then av- 
erage and normalise by the time interval 



k ± (t,At) = ± t (K±(t,At)). 



(2) 



Here and throughout, averages (...) run over the stochas- 
tic dynamics of the system and over a distribution of 
initial configurations where particle positions are chosen 
independently at random. In the limit At — > 0, then 
the fc ± (i,Ai) converge to the rates for bond-breaking 
(+) and bond-making (— ) events, which we denote by 
fc ± (t). The flux, f(t) = k+(t) - k~(t) is the net rate of 
bond-making and the traffic is the total rate of all events 
r(t) = k + (t) + k~(t). (In contrast to i23] . we focus here 
on rates for making and breaking bonds and not on the 
total numbers of bonds made or broken.) 
We also define a flux-traffic ratio 



Q(t) 



r(t) 



(3) 



which provides a dimensionless measure of the instanta- 
neous reversibility of a system, consistent with the qual- 
itative description of reversibility due to Whitesides [18J. 
The inverse of the flux-traffic ratio 1/Q(t) is the number 
of energy changing events per net bonding event. In a 
system at equilibrium f(t) = and so Q(t) — also. 
For a system quenched to T — 0, our MC dynamics does 
not permit bonds to be broken so particles never increase 
their energy: hence k~(t) = and Q(t) = 1. In systems 
that have been quenched, we expect to see a value be- 
tween these two limits, < Q(t) < 1: the smaller the 
flux-traffic ratio the 'more reversible' the system, while 
large flux-traffic ratios permit more rapid assembly. 
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FIG. 3: Plots of the flux f(t) and traffic r(t), and their ratio Q(t). The flux varies little across the temperature range 
considered while the traffic depends much more strongly on the temperature. In the kinetic trapping regime T/eh < 0.2, 
the ratio Q remains relatively large throughout the trajectories, while the good assembly and high-temperature regimes are 
associated with reversible evolution and small values of Q. 



In Fig. [3] we plot the flux f(t), traffic r(t), and their 
ratio Q(t) at a range of bond strengths, 2 < eb/T < 10. 
In all cases, the flux decreases towards zero as the system 
evolves towards equilibrium. On this logarithmic scale, 
the temperature-dependence of the flux appears quite 
weak, although the difference between different bond 
strengths may be up to an order of magnitude. Given 
that the 'integrated flux' F(t) = f*dt'f(t') p3](Figj3k 
inset) is always less than the maximal number of possi- 
ble bonds (4 in this case), it is clear that f(t) < i -1 at 
long times. 

In contrast to the flux, the traffic shows a large varia- 
tion with bond strength. Weaker bonds are more easily 
broken and result in more traffic. Combining flux and 
traffic, the ratio Q(t) is less than 0.1 throughout the 
good assembly regime, and decreases approximately as 
t^ 1 while good assembly is taking place. The regime of 
kinetic trapping is characterised by larger values of Q(t) 
and weaker time dependence. Similar results were found 
in Ref. [23 , where a similar ratio denoted by M was used 
to compare total numbers of bond-breaking and bond- 
making events. Compared to M, the ratio Q(t) depends 
only on the state of the system at time t and not on its 
history: this will be useful in making contact with other 
measures of reversibility to be discussed below. The ratio 
1/Q(t) is in the range 10 — 1000 during optimal assem- 
bly, indicating (as in [23]) that the system makes many 
steps forwards and backwards before it achieves a single 
step of net progress towards the assembled state. In this 
sense, good assembly requires significant reversibility, as 
argued above. 



B. 'Flux relation' between correlation and response 
functions 

In designing and controlling self-assembly processes, it 
would be useful if reversibility could be measured and 



controlled, to avoid kinetic trapping effects. While flux 
and traffic observables are not readily measured except 
in simple simulation models, we now show how correla- 
tion and response functions can be used to reveal simi- 
lar information (see also [H)l 122*]). These functions can 
be measured without the requirement to identify specific 
bonding and unbonding events; in some cases correlations 
and responses can also be calculated experimentally [34- 

cei. 

At equilibrium, fluctuation dissipation theorems 
(FDTs) [37] allow the response of a system to an exter- 
nal perturbation to be calculated from correlation func- 
tions measured in the absence of the perturbation. Away 
from equilibrium, such comparisons may be used for a 
classification of aging phenomena [2S] and may also be 
useful for estimating the degree of (ir) reversibility in self- 
assembling systems [TUI |2"2"] . 

We first consider the general case for measurement of 
a response function in a system with MC dynamics. Fig. 
[4] illustrates our notation and the procedure used. The 
system is initialised at t = and evolves up to time 
w, when a perturbing field of strength h is switched on. 
A single MC move is attempted with the field in place, 
after which the field is switched off. The system evolves 
with unperturbed dynamics up to a final time t at which 
an observable A is measured. In general, the response 
function depends on the two times t and w and gives the 
change in the average value of A, in response to the small 
field h. As indicated in Fig. [4] we use Greek letters to 
represent configurations of the system. (Recall that 5t is 
the time associated with a single MC move.) 

The MC scheme may be specified through the proba- 
bilities P°(i/ -s— (j,) for transitions from configuration [i 
to v in a single attempted move. However, it is con- 
venient to work with transition rates W°{v <— fi) = 
(5t)~ 1 P°(i' <— jtt). We also define p w {n) as the proba- 
bility that the system is in the specific configuration fi 
at time w, and the propagator G^{y <— fx) as the prob- 
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FIG. 4: Procedure for measuring the impulse response. The 
system is initialised in a random configuration. After a wait- 
ing time, w, a perturbation is applied for a single MC move 
during which the system attempts a move from configuration 
/i to v. The field is then switched off and the system allowed 
to evolve until time t, at which point its configuration is de- 
noted by 7. The typical response of the system is indicated, 
together with snapshots of showing how particles might move 
through the system. (The cluster that moves in the step while 
the perturbation is applied is highlighted in red.) 



ability that a system initially in state p will evolve to 
state v over a specified time period t. The superscripts 
indicate that no perturbation is being applied to the 
system: we use a superscript h if the field is applied. [It 
follows that p w {p) = XL^tuCA* K )Po( K ) where Pq(k) 
is the probability of initial condition k in the dynamical 
simulations.] 

With these definitions, the average value of the observ- 
able A at time t is 

(A(t)) = £ A( 7 )G U ( 7 <- v)P h (y <- p)p w {p), (4) 



where A (7) is the value of A in configuration 7 and we 
introduce u = t — w — St for compactness of notation. 
(Since we use Greek letters to represent configurations, 
there should be no confusion between ^.(7) and (A(t)), 
the former being a property of configuration 7 and the 
latter a time-dependent average.) The definition of the 
(impulse) response is 



R(t,w) 



T d(A(t)) 



(5) 



St dh ' 

where the derivative is evaluated at h — 0. Hence, 

R(t,w) = T^A( 7 )G°( 7 <- v)-^W h (v <- p) Pw (p). 

(6) 

We assume that the system obeys detailed balance 
with respect to an energy function E = E° — hB where 
E is the energy of the unperturbed system and B is the 



conjugate observable to the field h. It is convenient to 
define a connected correlation function 



C(t,w) = (SA(t)SB(w)} 



(7) 



(here and throughout we use the notation SO = O — (O) 
for all observables O). We also define 



S(t,w) 



_d_ 

dw 



C(t,w). 



(8) 



[The MC dynamics evolves in discrete time steps: the 
interpretation of the time-derivative is discussed in Ap- 
pendix lAl] At equilibrium, one has the FDT 



R cq (t,w) = S cq {t,w). 
Out of equilibrium, we define 

A(t,w) = S(t,w) - R(t,w) 



(9) 



(10) 



as a deviation between the response of the actual system 
and the response of an equilibrated system with the same 
correlation functions. 

Various out-of-equilibrium expressions for response 
functions may be derived: see for example 1551 159"] 
where different representations are defined and compared 
with each other. In Appendix A 1 we give a proof of ^ 
and we derive the general relation 

A(t,w) = J>( 7 ) - (A(t))]Guh <- u)F(u,M(p,n) 

(11) 



where 



and 



(12) 



J°G«, v) = W°(u «- p) Pw (p) -W°(p^- p)p w {v) (13) 

is the probability current between configurations p and v 
at time w. We will see below that the probability current 
J° plays a central role in quantifying irreversibility and 
deviations from equilibrium. 



The formula (11) is general for discrete time Markov 
processes obeying detailed balance (the generalisation to 
Markov jump processes in continuous time is straight- 
forward) . We emphasise that (11) is a representation of 
mj, valid whenever detailed balance holds. As such, it 
is mathematically equivalent to several other representa- 
tions that have been derived elsewhere: for example, the 
'asymmetry' function discussed by Lippiello et al. |39]I40] 
plays a role similar to A(i, w) in the case of Langevin pro- 
cesses. The purpose of (111 is to make contact between 



reversibility and deviations from FDT in self-assembling 
systems, as we now discuss. Further details of the rela- 
tion between (11) and previous analyses 1531 155HT21 are 
discussed in Sec. 



mG\ below. 
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FIG. 5: Data for SAt(t,w), t?At(t, w) and A A t(t,w) at three representative temperatures. The data are obtained at fixed 
w = f 3 MCS, while the time t is varied. The deviation between correlation S(t, w) and response R(t, w) is small in both good 
assembly and poor assembly regimes, but significant in the kinetically frustrated regime. We emphasise that the correlation 
and response functions are associated with perturbations to particle bond strengths: that is, A = B = n v in the definitions of 
S and R. 



C. Energy correlation and response functions in 
the lattice gas 

We now turn to response functions in the lattice gas. 
We consider how a single particle in the lattice gas re- 
sponds if the strength of its bonds are increased. To this 
end, we write the energy in the presence of perturbing 
fields h p as E h = J2 P ih eh ~ h p ]n p where the sum runs over 
all particles, as in (fil). We measure the response of (n p ) 
to the field h p , by taking the observables A and B of the 
previous section both equal to the number of bonds n p for 

a specific particle p. Thus R(t,w) = ^ d ^lf l ■ Given the 

energy function E h , there is still considerable freedom to 
choose the MC rates W h while preserving detailed bal- 
ance. 

In [551 [551 HO], it was mathematically convenient to 
take W h (y <- (j,) = W°(u <- ^) e h[B(p.)-B(u)\/2 _ Here 

our dynamics we are motivated by the central assump- 
tions (i) and (ii) of Sec. II A| which indicate that rates 
for bond-making and cluster diffusion should depend only 
weakly on perturbations to the particle bond strengths. 
We therefore include all /i-dependence in the probabil- 
ity of accepting an MC move, as follows. If the 
change in the unperturbed energy for an MC move is 
AE° and the change in the perturbation is AV^ = 
~J2 P h P An P then for AE ° / Owe take take P£ = 
min(l,e( A - 1 ) Af!0 / T - AV '/ T ) while for AE° = we take 
P% = 2a/ (1 + e Ay / T ). It is easily verified that this 
choice is compatible with detailed balance and reduces 
to the unperturbed probabilities P a if Ay = 0. Further, 
coupling the field h p only to the acceptance probability 
ensures that the perturbation affects the rates for bond- 
breaking and bond-making but does not affect the rates 
for diffusion of whole clusters. 



We use a straightforward generalisation of the 'no-field' 
method |4"Tl |4"2"] to allow efficient measurement of the re- 
sponse: see Appendix |A 2| for details. To attain good 
statistics, we consider responses in which the perturbing 
field h acts not just for one MC move but for a time inter- 
val At = 10 MC sweeps. Since we are working at leading 
order in h, the response to such a perturbation is simply 



St 



(At/St)-1 



R At (t,w) = — R(t,w + jSt). (14) 



3=0 



The relevant correlation function in this case is 
C(t,w) = (6n p (t)6n p (w)) where Sn p (t) = n p (t) — (n p (t)) 
as usual. It is convenient to define normalised correlation 
and response functions 

S At (t, w) = [C(t, w + At)- C(t, w)] (15) 



R A t(t,w) 



R A t{t,w) 



(16) 



with j\f At {w) = C(w + At,w + At) - C(w + At,w) so 
that S A t(w + At,w) = 1. We also define 



A At(t, w) = S At (t, w) - R At (t, w). 



(17) 



In all cases, the subscript At indicates that these func- 
tions generalise the R(t,w), S(t,w) and A(t, w) of the 
previous section, while the tilde indicates normalisation. 
For small At, the dependence of these functions on At is 
weak, but numerical calculation of the response is easier 
for larger At. 

In Fig. [5] we show numerical results for correlation and 
response functions. At T/e^ — 0.5, clusters of parti- 
cles are growing in the system, which is far from global 
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equilibrium. However, as found in |19l 122] for several 
other self-assembling systems, the deviations from FDT 
are small because the time-evolution is nearly reversible 
(recall Fig. [3|. On the other hand, at T/eb = 0.15, 
the particles are aggregating in disordered clusters and 
the time evolution is far from reversible with unbonding 
events being rare. 

The potential utility of this result was discussed 
in [TU1 122] : it means that straightforward measurements 
on short time scales can be used to predict long-time 
assembly yield, by exploiting links between correlation- 
response measurements and reversibility. In what fol- 
lows, we explore in more detail how the deviation func- 
tion A couples to microscopic irreversibility during as- 
sembly. 



D. 



Interpretation of A(w + St, w) as a flux 



In this section, we concentrate on quantities that can 
be written in the form 



Z(t) 



(18) 



where J° is the probability current, defined in ( 13 1 



In contrast to straightforward one-time observables like 
(A(t)) — ^^A^Ptin), we wm show that observables 
like Z(t) are currents, time-derivatives and fluxes: they 
measure deviations from instantaneously reversible be- 
haviour, measured at time t. 

To this end, we generalise the analysis of "kinks" in 
Sec. Ill A by defining (K "^(t, At)) as the average num- 
ber of MC transitions from state /i to state v between 
times t and t + At. The associated kink rate is k v ^{t) = 
(St)- 1 (K"^(t,St)). It follows from the definition of the 
fes that 



k u »(t)~k^(t) = J t V,/x) 



(19) 



Thus, the probability current J t ° gives the difference 
between the probabilities of forward and reverse MC 
transitions between states /z and v, which clarifies its 
connection to irreversibility. At equilibrium, one has 
J t °(i/, n) = for all states /i and v. thus, J t ° measures de- 
viations from equilibrium. However, measuring (or even 
representing) J t ° is a near-impossible task in a system 
where the number of possible states /i is exponentially 
large in the system size. 

Instead, one makes a specific choice of the z(y, fi), in 
which case ( 18 1 shows that Z(t) is a projection of the 
current J t ° onto the observables of the system. For ex- 
ample, if z(y,(j,) = A{v) then Z(t) = ^(A(t)) is a time 
derivative that clearly vanishes at equilibrium (or in any 
steady state). Other choices for z(i/,fi) become relevant 
in out-of-equilibrium settings. For example, the flux f(t) 
in Sec. Ill A is obtained by setting z{y, /i) = 1 if the tran- 
sition ji — >• v involves an increase in the number of bonds 
of particle p, with z(y, /a) = for all other /z, v. (The 



A At (t + At,t) a) 
l 



0.75 




0.25 



1500 



A At (i + AM) b) 
1 



0.75 




0.25 



Q(t) 



FIG. 6: a) The deviation from FDT measured immediately 
after perturbation, A& t (t + At, t), At = 10. At temperatures 
above T* , A&t(t + At,t) decreases quickly with t, while at 
low temperatures it remains significant throughout the early 
times. In the construction of A(t,w) we take the observables 
to be A = B — - Tip cLS discussed in the main text, (b) Para- 
metric plot of Aa*(£ + At, t) against the flux-traffic ratio Q(t), 
for < t < 1500. The relationship between these two quanti- 
ties depends very weakly on temperature, indicating the close 
relation between them. 



structure of (18 1 together with this choice of z ensures 



that f(t) acquires negative contributions from MC tran- 
sitions where particle p experiences a decrease in its num- 
ber of bonds, as required.) Similarly the particle current 
in exclusion processes is obtained by taking z(v,p) = 1 
for transitions /i — > v where a particle hops to the right, 
and zero for all other pairs. Clearly, these observables 
monitor the breaking of time-reversal symmetry, as is 
the case for all observables Z(t). 

To relate deviations from FDT to irreversible events, 
we first consider the case t = w + St, so that the per- 
turbation h is applied for just one MC move and then 
the response is measured immediately. In this case (111 
reduces to 

A(w + 5t,w)=^2zA(y,ii)J°{u,fi) (20) 

/A// 
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as in ( 18 ), with 



za(", M) = [My) ~ (A(w + St))]F{y, /x). (21) 



The function F[y, fi) was defined in (JT2J) and measures 
the effect of the perturbation h on the transition rate 
from v to p. The factor A(v) — (A(w + St)) indicates 
whether configuration v has a high or a low value of the 
observable A, compared to the average of A at the mea- 
surement time w + St. 

For the case considered in Sec. Ill C and Fig. [5j where 
the perturbation is coupled to the energy of particle p, 
one has A(v) — n p (v) and 

F(u, M ) = [n p (u) - n p (/i)]e(E°(/i) - £») (22) 

where Q(x) is the step function, defined such that O(0) = 
2 • The step function appears because if a proposed MC 
move results in a decrease in the bare energy, the ac- 
ceptance probability does not depend on the perturbing 




0.4 0.6 

C(t,w) 



FIG. 7: Integrated correlation-response plot showing the 
response x(t,w) and correlation C(t,w). The correlation and 
response functions are associated with observables A = B = 
as in previous figures. The data is shown for a fixed 



field h p . Clearly, F{y, /i) is finite only if particle p changes t = 10 4 MCS and and < w < t. The x symbol indicates 



its energy between states \i and v. the same is true for 

Thus, (20 1 indicates that the deviation A(t,w) mea- 
sured at t = w + St reflects the imbalance between 
rates of bond-making and bond-breaking for particle p. 
In this respect it is similar to the flux f(w): however, 
the weights given to different bond-making and bond- 
breaking processes differ between A(w + St,w) and f(w), 
due to the different forms of z(u, p) in the two cases. 
Thus, while these quantities reflect similar physics, the 
details of their behaviour is different. To see the rela- 
tionship between the immediate response A(w,w + St 
and the flux f(t) from Sec. Ill A we present Fig 
which shows that the two quantities have similar time 
and temperature-dependence and therefore reveal similar 
information about the (ir)reversibility of the dynamical 
bonding and unbonding processes. 



E. Time intervals for reversibility 

We now consider the deviation A(t, w) for t > w + St. 
From (11 1, we see that A(t, w) can be written in the form 



( 18 1 if we take 



where 



z(v,(x) = [V A (v)-(A(t))]F(v,»). (23) 



P»=5>( 7 )G°(7^) (24) 



the points where w = 10 : impulse responses associated with 
this time are shown in Fig. [5] 



since the system forgets its memory of the configuration 
v and "regresses back to the mean". (That is, the propen- 
sities for different fixed initial conditions all converge to 
the same average value at long times.) Out of equilib- 
rium, this may not be the case: if a configuration v has 
a more ordered structure than another configuration v' 
then these configurations may have different propensities 
for order even at much later times. Loosely, a system in 
v has a 'head start' along the route to assembly, and the 
system does not forget the memory of this head start as 
assembly proceeds. 

Fig. [5]shows that for the correlation and response func- 
tions considered here, the deviation A(t,w) decays only 
very slowly with time t. The f-dependence of this func- 
tion comes through a weighted sum of propensities: it 
is therefore clear that these propensities do not regress 
quickly to the mean. In contrast, one may write the cor- 
relation as 

S(t,w) = Ypuiy) - (A(t))]\nM - n p (p)} 

xW°(is^ ii)pM (25) 



is the propensity of observable A for configuration v 
after a time u. That is, T^iv) is the value of A that is 
obtained by an average over the dynamics of the system, 
for a fixed initial condition v and a fixed time u. (Recall 
u = t — w — St.) 



Comparing (|23[) with (21 1, the difference is the replace- 



ment of A(v) by the propensity. At equilibrium, one ex- 
pects V^(y) — (A(t)) to decay towards zero as u increases, 



whose i-dependence also comes from the propensity. 
Fig. [5] shows that this function does decay quite quickly 
with t, although it does not reach zero on the time scales 
considered here. 

Our conclusion is the following. For configurations v 
that are typical at time w, the propensity V^iy) has a 
'fast' contribution that decays quickly with time u, as 
well as a 'slow' contribution that depends weakly on u 
and reflects the 'head start' of v along the route to the 
assembled product. The fast contribution reflects rapid 
bond-making and bond-breaking events that do not lead 
directly to assembly while the 'slow' contribution is a 



9 



non-equilibrium effect that measures assembly progress 
and also dominates the deviation A(t,w) that we have 
defined here. By contrast S(t,w) picks up contributions 
from both slow and fast contributions. The utility of 
the FDR is that for suitable w, it couples to the irre- 
versible bond-making behaviour that is most relevant for 
self-assembly. 

To illustrate this balance of fast and slow degrees of 
freedom, and also to make contact with previous stud- 
ies [T9] 122] [26] of FDRs, we define an integrated (and 
normalised) response 



X{t,w) 



1 



C(t,t) 



dw' R(t, w') 



and we also normalise the correlation as 



C(t,w) 



C(t,w) 
C(M7 



(26) 



(27) 



It is conventional to display these correlations func- 
tions in a parametric form [26;. Specifically, making a 
parametric plot of x(i, w) against C'(t, u>), for fixed t and 
varying w, the gradient is —X(t,w). It has been empha- 
sised several times pHl - HS] that plotting data for fixed w 
and varying t is not in general equivalent to this proce- 
dure and may give misleading results, especially if C(t, t) 
has significant time-dependence, as it does in these sys- 
tems. Plotting data at fixed t as in Fig. [7] it is ap- 
parent that the high temperature systems are the most 
reversible, with X(t,w) « 1. To understand the con- 
nection with Fig. [5j note that the time interval t — w 
increases from right to left in Fig. [7j when t — w is large 
then the deviation A(t,w) has a larger fractional contri- 
bution to S(t,w), so that the curves are steepest (most 
reversible) when t — w is small and least steep (least re- 
versible) when t — w is large. We find that parametric 
plots depend weakly on the fixed values of t used but 
we do not show this data, for brevity: see for example 
Ref. [22] for an analysis of i-dependence in a system of 
crystallising particles. 

Physically, we attribute the w-dependence of X(t,w) 
to the fact that fast degrees of freedom (like dimer for- 
mation and breakage in the vapour phase) quickly relax 
to a quasiequilibrium state [27] where probability cur- 
rents J° associated with this motion are small. It is 
only much slower degrees of freedom (like the gradual 
growth/assembly of large clusters) for which probability 
currents remain significant on long time scales, and lead 
to long-time contributions to A(t,w). 

Compared to the simple flux-traffic ratio Q(t) consid- 
ered in Sec. Ill A the normalised deviation 1 — X(t, w) = 



A(t,w)/S(t,w) has the same effect of comparing a gen- 
eralised flux A(t, w) with a normalisation S(t, w) that re- 
flects the traffic in the system. However, the effect of the 
time difference t — w is that while A(t, w) decays slowly 
with time t then S(t, w) decays quite quickly. Thus, as 
t increases, the deviation 1 — X(t, w) becomes increas- 
ingly sensitive to the deviations from reversibility, since 




0.4 0.6 

C{t,w) 

FIG. 8: Integrated correlation-response plot a site dependent 
perturbation. That is, we take the observables A = B = pi as 
described Sec. |IHF] and in contrast to previous figures where 
A = B = n p . The data is shown for fixed t — 10 4 with 
x symbols indicating w — 10 3 , as in Fig. [7] The relaxation of 
the observable pi is diffusive and depends weakly on the bond 
strength, so the dependence on T/tb is much weaker than in 
Fig. [7] and the deviations from FDT behaviour are smaller. 



the fluctuations of fast (reversible) degrees of freedom 
regress back to the mean while the slow (irreversible) ones 
remain significant. For this reason, the FDR X(t,w) is 
a more sensitive measure of irreversibility than the flux- 
traffic ratio discussed in Sec. IIII A"l 



F. Responses to different perturbations 

We remark that the observables used in measuring cor- 
relation and response functions may strongly affect the 
results. Most relevant is the extent to which the response 
function couples to the slow, irreversible degrees of free- 
dom and the fast, reversible ones. To illustrate this, we 
compare the results of the previous section with a differ- 
ent response function. We take the observables A and B 
of Sec. |III B| both equal to the occupancy pi = 0, 1 of a 
specific site i on the lattice. The effect of the perturba- 
tion onthe MC dynamics is the same as that described in 
Sec. IIIC except that the perturbed energy E h = E° + V 
with V = J^ihiPi, where the sum runs over sites of the 
lattice. 

Unlike the bond numbers n p which change only when 
bonds are made and broken, the site occupancies pi also 
change as clusters diffuse through the system. When 
bonds are strong, the site correlation functions decay 
faster than the bond correlation functions considered 
above. Therefore, since the site response is coupling to 
fast degrees of freedom, we expect to see 'more reversible' 
behaviour: this is borne out by Fig. [8] where the FDRs 
are much closer to the equilibrium value (unity) than the 
FDRs shown in Fig. [7] This is consistent with recent re- 
sults of Russo and Sciortino |46| who observed a similar 
effect, and with earlier studies of observable-dependence 
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of FDRs [44 451137] (see however gSj in which some fast 
degrees of freedom do not relax to quasiequilibrium) . 

Similar correlation and response measurements in an 
Ising lattice gas were also considered by Krzakala |49| . 
but the dynamical scheme used in that work means that 
clusters of particles do not diffuse: in that case the 
site degrees of freedom pi are not able to reach quasi- 
equilibrium and one does observe significant deviations 
from FDT for this observable. 

We end this section by briefly considering these corre- 
lation and response functions for large t. In this limit, 
standard arguments [25] indicate that both correlation 
and response functions have contributions from well- 
separated fast and slow sectors, with the contributions 
from the fast sector obeying FDT. In coarsening systems 
like this one then there is no response in the slow sec- 
tor. Thus, for large enough times, the parametric plot 
has a segment with X = 1 for large C and small \'i 
for smaller C then X = so the parametric plot has a 
plateau where the value of \{t, w) depends only on equi- 
librium properties of the system. However, we emphasise 
that the approach to this large-t limit is very slow and 
applies only after all clusters in the system are annealed 
into compact clusters. This limit is therefore irrelevant 
for analysing the kinetic trapping phenomena on which 
we focus in this article. 



G. Relation to previous analyses of 
non-equilibrium response 

We have emphasised throughout this article that our 
main purpose is to use correlation and response mea- 
surements to measure reversibility in assembling systems. 
Many other studies have considered such measurements 
in a variety of other contexts - here we make connections 
between our methodology and some other results from 
the literature. 

One recent area of interest has been the use of 'no- 
field' measurements of response functions. Here, the aim 
is to develop formulae for the response that can be eval- 
uated in Monte Carlo simulations without the introduc- 
tion of any direct perturbation |39f[4"2"l 14"5], In fact, we 
use such a straightforward generalisation of the method 
of |41| to measure responses in this paper, as discussed 
in Appendix |A 2| 

Our analysis in Sec. |III B| is concerned not so much 
with measurement of the response, but with the phys- 
ical interpretation of deviations from FDT. Our result 
is therefore in the same spirit as [3U |2"1H21)1 I53"l 1551 14"0] . 
In particular, Baiesi et al [33J relate the FDT to 'flux' 
and 'traffic' observables that separate reversible and ir- 
reversible behaviour, but they use a different approach 
to ours. Full details are given in Appendix. |A 3| but in 
our notation, their central result is 



where T(w) measures the /i-dependence of the amount 
of dynamical activity (traffic) between times w and w + 
St. The second term is therefore associated with traffic 
while the first term is a correlation between A{t) and 
the 'excess entropy production' at time w: it is therefore 
related to a flux. However, neither of the correlation 
functions in (|28| may be written in the form of ( 18 1 so 



they do not vanish at equilibrium. The key point is that 



in equilibrated (reversible) systems, the terms in (28 1 are 
both non-zero and equal to each other. So while (28 1 



does involve a separation of terms symmetric and anti- 
symmetric under time reversal, the resulting correlation 
functions are not fluxes in the form of (18) and do not 



provide the same direct measurement of the deviation 
from reversibility given by A(t,w) in (11 1. 



In another relation to irreversible behaviour, the anal- 
ysis of Refs. |391 14TJ] identifies an 'asymmetry' measure- 
ment which in our notation is 

A(t,w) = R(t,w) - -[S(t,w)+S(w,t)] (29) 

It follows from the fluctuation dissipation theory that this 
quantity vanishes in a system with time-reversal symme- 
try. Comparison with (10 1 also illustrates the similar- 



ity between A(t,w) and the deviation A(t,w). However, 
A(t, w) differs from A(t, w) since the deviation A(t, w) is 
a projection of J° and hence vanishes if the system is be- 
having reversibly at time w, regardless of what happens 
at later times. On the other hand, A(t, w) vanishes only 
if all trajectories between times w and t are equiproba- 
ble with their time-reversed counterparts. That is, time- 
reversal symmetry must hold throughout the interval be- 
tween w and t, and not just at time w. 



IV. OUTLOOK 

We have analysed the extent to which ideas of re- 
versible bond-making are applicable to a phase separat- 
ing lattice gas, which we used as a simplified model for 
a self-assembling system. The data of Figs, [I] [3] and [JJ 
further support the idea that this system is relevant for 
studying self-assembly, since similar results have been 
presented for more detailed models of self-assembling sys- 
tems |9j[19l|20l|22l[23]. 



The flux/ traffic measurements of Sec. Ill A provide an 
intuitive picture of reversibility, while the picture based 
on the correlation-response formalism is more subtle. 
However, the central result (111 demonstrates an explicit 



R{t,w) 



1 



d_ 
dw 



(A(t)B(w)) + (A(t)T(w)} 



(28) 



link between measurements of response functions and ir- 
reversibility of bonding. Also, Figs. [5] and [JJ do indicate 
that information about both short-time reversible and 
long-time irreversible behaviour can be obtained by con- 
sidering the behaviour of correlation and response func- 
tions. 

It has been argued that correlation-response measure- 
ments on short time scales might be used to predict long- 
term assembly yield, and even to control assembly pro- 
cesses. The work presented here places this objective 
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on a firmer theoretical footing, and highlights the im- 
portance of using the right observable when probing irre- 
versibility (compare Figs. [7] and [8]); the interplay between 
fast 'quasiequilibrated' degrees of freedom and slow 'ir- 
reversible' bonding [27 also highlights the importance of 
choosing the right time scale for measuring these func- 
tions and predicting long-time assembly quality. 
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Appendix A: Derivations of representations for 
response functions 

1. Deviation function A(t, w) 

In this Section, we derive equation ( |11[ ) which gives 
the deviation between the response R(t/w) and the cor- 
relation function S(t,w). Starting from ([6|, conservation 
of probability for the MC transition probabilities implies 
that P h (n <— fi) = 1 — X^i/fAi) P h { v <— A 1 ), with a similar 
relation for the rates W h . Using this in ^ gives 



R(t, w)=Tj2 Ml)Gu(l <- v)^Jt{^ v) (Al) 



with 



v) = W h {v <- ripvfji) - W h (fi «- v) Pw {v) (A2) 

which has the structure of a flux in probability from con- 
figuration fi to configuration v. This is related [38 to an 
expression for the response originally due to Agarwal [5U] . 

We now assume that the system obeys detailed bal- 
ance: 

Whfr <_ v ) e -&(?)/T = W h {v ^ fl)e -E h M/T (A3) 

where E h ([i) = E°(fi) — hB(^i) as defined in Sec 
It 

d 



Hence, one has ■^lnW h (i' <— //) 



IIIB 



M \nW h ((i 4- v) with AB(u,fj,) = B{v) - B(fj,), so that 



+ T.f w { V ^)-^\nW h ( i x<-u). (A4) 



Inserting ( A4| into ( Al ) gives two contributions to the 
response. The first contribution is ir(A(t)[B(w + St) — 
B(w)]}. In the notation of the main text, this is equal 



to S(t, w) + (A(t)) (B(w)) where the time derivative is 
interpreted as a change between w and w + St, normalised 
by the time St for a single MC move. (Recall that we took 
time to be discrete in these MC models.) 

At equilibrium, J° — and -J^(B(w)) — so the only 
contribution to the response is S{t, w) and the FDT ((9} 
holds. The non-equilibrium contributions to R(t, w) can 
be collected together to obtain (111 if one notes (i) that 
that ^ 7 G° (7 <— v) = 1 which follows from conservation 
of probability and the definition of G° and (ii) that 



d_ 
dw 



[B(w)), 



(A5) 



which follows from the detailed balance property of the 
W h . 



2. 'No-field' formula for the response 

Another useful representation of the response function 
expresses it in terms of directly observable quantities in 
unperturbed simulations |41ll^2"] (this is the 'path weight' 
representation of |38|i. 

The original 'no-field' method requires an extension for 
our system since the cluster MC algorithm we use means 
that the same transition v 4— fj, may take place via several 
possible computational routes. (For example, two MC 
moves may start from different seed particles but result 
in the same cluster being moved in the same direction). 

We use P(y < — fj,) to represent the probability of an MC 
move from fi to v by some route C. For v ^ /x, the route 
C is the choice of seed particle and sequence of steps 
by which the moving cluster is generated. If v — fi the 
route C may involve the proposal of a move that is then 
rejected, or the proposal of a move where the cluster size 
exceeds n max , or would result in multiple-occupancy of a 
site. 

Starting from ^ , one writes 

R(t,w) = T J2 AW&jru-vyn&CiJL) 

xP°(i/^^0*). (A6) 

where K{v,C,n) = j- t -§g log W h (y 4^ /1). This result 
may be written MS i\ simple expectation value 



R(t,w) = T(A(t)K(w)) 



(A7) 



where the key point is that the observable lZ(v, C, fi) may 
be evaluated for any attempted MC move so that the cor- 



relation function in ( A7 1 may be calculated directly from 



an MC simulation with h — 0. [This is in contrast to the 
current (p, /x) which depends on the whole ensemble 
of evolving systems, so that the formula ( Al I cannot be 
evaluated directly by MC simulation.] 

For the dynamical method described in this paper, 
the observable TZ(f,C,fi) is given by 7Z — [n p (v) — 
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n p (/i)]0(£' o (^) — E°(v)) if the move is accepted, as in 
the main text 0(0) — \. If the move is rejected ei- 
ther because the cluster size n exceeds n max or the move 
would have resulted in overlapping particles 1Z = 0. If 
however the move is legal but rejected when testing the 
move acceptance probability P a , then 1Z = —[n p (y) — 
n p (/i)]e(£ ( M ) - P») x P a /(1 - P a ). 



3. An alternative formula for the response 

In this section we clarify the connections between our 
results and those of Baiesi et al . 1331 . The central identity 
in the derivation of Appendix 



AB(i/, n)/T 
write 



Al 



is ^\nW(v <- M ) 



lnlE(/i -s— v). Rearranging, one may 



dh 



\yiW[v 



+ \iHW{n^v)W{v^n)] (A8) 



2dh 

and substitution into ^ leads to 

x [AB(v, n) + 7>, ti)]W(v <- mK(m) (A9) 



with T(i/, m) = ln[W(> <- v)W{v <- /x)]. Equ. (|A9| 
is the translation of the central result of Baiesi e£ a/. 733" 
into our notation. 

The key point here is that the two terms in square 



brackets in (A9l have different symmetry properties. By 



definition of the first term, AB(y,/t) = —AB(fi,i/). 
Time-reversal of the trajectory flips the roles of v and fi 
so that AP is 'anti-symmetric under time-reversal' and 
is related to the entropy current [33]. This means that 
the average of AP is zero in a time-reversal symmet- 
ric system. (Note however that AP itself is a fluctu- 
ating quantity and not zero in general: this is different 
from the probability current which is an ensemble- 
averaged property without fluctuations and is strictly 
equal to zero if time-reversal holds.) On the other hand 
T(y, n) — 7"(/Lt, v) so T is 'symmetric under time-reversal' 
and measures dynamical activity, or traffic. 



Despite the links between Equ. (A9) and time-reversal 
symmetries, neither of the terms in (A9l may be written 
as a flux in the form of (18 1. Noting that S(t,w) is a 



connected correlation function, one may write 



R(t,w) 



S(t,w) + (A(t)T(w)) 



(A(t))£(B( w) ) 



(A10) 

where the final term vanishes in the non-equilibrium 
steady states considered in ]33 while T(w) is given by 
(St)~ 1 7~(is, At), evaluated with v and \i being the configu- 
rations at times w + St and w respectively. The key point 



is that the first two terms on the right hand side of ( A10 1 



are equal to each other and both are non-zero at equilib- 
rium. Hence neither of them is directly identifiable as a 
flux that is sensitive only to deviations from reversibility. 
Combining the first and third terms on the right hand 
side of (A10) yields (28 1 of the main text. 
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